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The arrival time distribution of cosmic ray events is well suited to extract information regarding sky anisotropics. 
For an experiment with nearly constant exposure, the frequency resolution one can achieve is given by the in- 
verse of the time T during which the data was recorded. For T larger than one calendar year the resolution 
becomes sufficient to resolve the sidereal and diurnal frequencies. Using a Fourier expansion on a modified 
time parameter, we show in this note that one can accurately extract sidereal modulations without knowledge of 
the experimental coverage. This procedure also gives the full frequency pattern of the event sample under stud- 
ies which contains important information about possible systematics entering in the sidereal analysis. We also 
show how this method allows to correct for those systematics. Finally, we show that a two-dimensional anal- 
ysis, in the form of the spherical harmonic (Y™) decomposition, can be performed under the same conditions 
for all m ^ 0. 



1. Introduction and Notations 

If anisotropics of the incident cosmic ray flux depend on right ascension, they induce on the observed event 
rate a modulation with a period equal to one sidereal day. If other sources of modulations (weather effects, 
detector acceptance, etc), be them known or unknown to the experimenter, exhibit some periodicities, they 
have no reasons to be related to the sidereal period. In general, the extraction of large scale anisotropy rely on 
the difficult task of computing the detector coverage map, i.e. the detector response to a perfectly uniform sky. 
This estimation is often delicate as unknown or ill-corrected effects may introduce systematics in this map, 
limiting the precision on the genuine anisotropics from a given data set. 

The Fourier transform over an exposure time T allows to separate frequencies with a resolution of the order 
of 1/T. Data acquired over a few years, even with variable detection conditions, can be analyzed using this 
transformation allowing for a clean separation between the sidereal frequency and the diurnal frequency, which 
differ by 1 calendar day over a calendar year. This separation can be obtained totally independently of the 
knowledge of the detector coverage map, hence avoiding the difficulty of its estimation. At higher frequencies, 
when looking at multiples of the diurnal or sidereal frequency the resolution is even better since the separation 
between the n th harmonics is n times larger than between the fundamentals . 

The Fourier transformation is an ideal tool to extract the frequency patterns from a set of event times. We 
define the (normalized) Fourier coefficients a(f) of set of times {i,} at frequency / (cycle uu = 2nf) as : 

a(f) = c(f)+is(f) (1) 

with 
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c(/) = w^2 Wi cos ^ 

i 
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S{f)= W i Sm ( 3 ) 
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with W — J2i w i an d where the sums run over all the TV events whose time distribution is studied. The event 
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weights Wi may incorporate known modulations of the data that one wants to remove from the expansion. 
When no correction is desired lUj = 1 Vi, and W = The amplitude at frequency / is then given by 
K/)l = V s2 (f) + c2 (/) while the P hase is defined as <j>(f) = tan" 1 (s(f), c(/)). 

In general cosmic ray surface detectors have a large instantaneous aperture and can observe a sizeable fraction 
of the sky at any given time. For example, a surface array whose detector can see air showers with zenith up 
to 60 degrees sees instantly l/4th of the whole sky. To be able to observe modulations on smaller scale it is 
therefore necessary to correct the time entering in Eqn. 1 and 2 by the right ascension phase of the event with 
respect to the Local Sidereal Time (LST, the right ascension of the stars above the local meridian). We then 
modify ti entering in the equations above as: 

= U + (RAi - LST(U)) (4) 

where LST(U) is the local sidereal time corresponding to UTC time ti and RAi is the right ascension of the 
event. Note that LST(^i) — RAi. With such a definition the Fourier Expansion in Modified Time (FEMT) 
(with unit weights Wi = 1) is equivalent at the sidereal frequency to the Rayleighanalysis. Note that without the 
introduction of this modified time an eventual high order (small angular scale with respect to the instantaneous 
aperture) modulation in sidereal time would average almost to zero. On the other hand an eventual diurnal 
modulation will be slightly reduced when using this parameter. 

2. Properties 

In this section we show how the Fourier expansion allows to completely separate, at first order, sidereal time 
modulations from experimental modulations without the need to correct for the detector coverage map. 

A given event population P(d, £) with declination d is proportional to : 

P(d,£) = kF(LST(g),d) [ A(6,<p,t)6(£-t)6{d- d(6,ip)dtd cos 0d<p (5) 

Jt,$,<p 

i.e. 

P(d,0 = kF(LST(0,d) A'{d,t + nT sld ) = k'F{ra,d)A s {d,t) (6) 

n= — oo 

where k, k' are constants, F represents the sky flux in the direction (ra, d) and A (A 1 , A s ) is the detector 
acceptance corresponding to our event selection and T s id is the sidereal period. 

The Fourier transform, P of P is the convolution product of the Fourier transform F of F and A s of A s : 

P(d, f) = K J df'F(d, f)A s (d, f - f) (7) 

Because F represents the sky flux it is a perfectly periodic functions of time (or of Q and F is given by an 
harmonic series of the fundamental sidereal frequency = 1 /T S id : 



1 We will place ourselve in this case except in section 4 
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F(dJ) = Y / a m (d)S(f-mf sld ) (8) 

m 

with ao — 2 due to our definition of the Fourier transformation above. 

To simplify the notation we will also assume that A s is an harmonic series but of the diurnal frequency f d '■ 

i(d,/) = ^6 m (d)5(/-m/ d ) (9) 

m 

again with b^ — 2. A may contain other fundamental frequencies such as f y and also be non periodic^] The 
one important assumption we make is that f sid does not enter as a special frequency in A as might f d or f y . 

We can assume that the modulations introduced by a m and b m with m ^ are small (the acceptance is 
relatively stable in time and the true anisotropics, if any, are small). Then, at first order in a m and b m , P 
reduces to: 

P(d, f) ~ 2 + £ a m (d)S(f - mf sld ) + J2 b m (d)S(f - mf d ) (10) 

m > m > 

We observe a decoupling of the modes. The relative modulations at frequency / in each declination band d is 
directly given by the modulus of P. In particular for each sidereal harmonic m we have : 

\P{d.mf sid )\ = \P{d.m)\ ~ Vc2(mf sld ) + s 2 (m/ aid ) (11) 

Here c and s are given by eqn. 2 and 3 above but with the sum running over all events with declination d. 

Moreover, the knowledge of the evolution of P(d, m) as a function of the declination d allows for the recon- 
struction of the corresponding multi-polar coefficients under certain conditions; if F is expressed as: 

F ( 1 + Yl ^irnYr^/2 - d, ra) I (12) 

\ (Z,m)^(0,0) / 

the function a m (d), for a given value of m, gives J2i> m a imP'^ l {^ /2 — d), where the "P™ are the associated 
Legendre polynomials. If the expansion may be truncated to a certain order, we can deduce from these con- 
straints the decomposition of F over the Y{ n (m ^ 0), independently of the knowledge of A (but the precision 
is limited by the available statistics, of course). 

In practice experimental effects will of course partly blur this ideal picture, in particular : 



• The sum in eqn. 2 and 3 only runs over a finite amount of time T hence the Fourier transform is 
convoluted with the transform of the function describing the acceptance as a function of the time. With 
a constant acceptance it is a function (smx/x) with x = ujT; it is more complex when accounting for 
possible overall increase or short scale irregularities of the detector up time, but in any case, with T larger 
than a few years it has a width less than the difference between the sidereal and diurnal frequencies. We 
will see below that this resolution is sufficient to perfectly separate the two components. 



2 A detector only active during a period of time T is a square function and its Fourier transform would be the convolution of 1/ (2nfT) X 
sin (27r/T) with the harmonic development above. 
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• We have a finite number of events N, hence a Poisson noise in all frequencies below about N/T, this will 
induce an average statistical noise on our amplitude estimates of order ^/ttJN. 

• The diurnal modulation are not perfectly periodic with a 24 hour period. Moreover annual variation over 
a period of about a year (355.25 * fd) also exist. This may produce second order (the amplitude is the 
product of the annual and diurnal amplitude) side bands at fd + f y — f s %d an d fd — fy — fanti-sid- In 
practice the existence of such bands can be checked on the anti-sidereal (lower frequencies) side of the 
frequency distribution. Therefore, after a first order reconstruction of F and A s from eqn. 10, one can 
evaluate the contribution of the second order terms, and, if necessary, to re-evaluate F correcting for the 
principal modulations of A s , An example of such corrections is given in section |4] 

3. Noise statistics and upper limits 

In the case of a totally uniform sky and in the absence of local effect the time entering in equation 1-3 are 
uniform and the variable c(/) and s(f) converge, in virtue of the central limit theorem, toward Gaussian 
random variables with mean and variance 2/N, Due to the nature of c(/) and s(f) the convergence is very 
rapid and holds in practice for N of order a few tens. The random variable z(f) — \J c(/) 2 + s(/) 2 is then a 
Rayleigh random variableof parameter a = W2/N whose distribution is : 

f(z | y/2/N) = ^ e -^ 2 /4 (13) 

The mean of this distribution is •v/^ and the variance is ^jf 1 - The cumulative distribution function giving the 
probability for the noise to be above some power zq is given by : 

F(z | y/2/N) = 1 - e- Nz °' A (14) 

Hence, as in the standard Rayleigh analysis of cosmic ray data, for a background sample of size N the a CL 
upper limit on the noise power is given by p a = — 2-y/ln (1 — a) /N. 

In practice what is often needed is, given a measured power p, the probability to observe this measurement 
or less when a non random modulation of power po is present in the data set. When signal is present the 
distribution of z is no longer given by Eqn. [T3]but by the Rice distribution : 

R(z | VWN,Po) = ^e- N ^ + ri» 4 I (z Po N/2) (15) 
where Iq(x) is the modified Bessel function. 

The a confidence level upper limit p a (/) on a signal at frequency / when a power p has been observed can be 
defined as the solution of: 



r+oo 

/ R(z | ^2/N, Pa )dz = a 
Jp 



(16) 
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Data set 


N 


Py 


Pd 




& stat 


side band 


diurnal resolution noise 










(v/tt/AO 


(V(4-tt)/JV) 


(Py X p d )/2 


~ Pd/(nT) (T in year). 


A 


10 5 


19 


4 


0.6 


0.3 


0.5 


0.4 


B 


4 x 10 5 


19 


4 


0.3 


0.15 


0.4 


0.3 


C 


2 x 10° 


19 


4 


0.1 


0.07 


0.4 


0.1 



Table 1. Background parameters of the different data sets. In data set A the statistical noise dominates, in data set B 
statistical and systematical noise are similar, in data set C the systematics dominate. Note that the amplitude given for the 
annual (p y ) and diurnal (pd) modulation are the one reconstructed in modified time hence they are smaller than the original 
ones in UTC time (20% and 5% respectively). The diurnal resolution noise is given at the sidereal (and anti-sidereal) 
frequency. All numbers, except N, are in percent. 



which can be integrated numerically. 

When the systematical noise, such as the resolution noise or the side-band noise, becomes dominant over the 
statistical noise the distribution of the reconstructed power at the sidereal frequency is only a Rice distribution 
for a fixed signal phase with respect to the systematical noise. When the signal phase is random the power dis- 
tribution becomes a superposition or Rice distribution with different po parameters (depending on the relative 
phase between the systematical noise and the signarb. 



4. One dimensional examples 

To illustrate the method and the various properties depicted in the previous section we constructed random data 
sets according to 3 different situations. All data sets have an annual modulation of 20% with a fixed phase 
4> y = 1.0, a diurnal modulation of amplitude 5% and phase (p^ = 0.0 and eventually either a dipole modulation 
of amplitude 1% declination 0° and phase ir/2 or a right ascension modulation of amplitude 1.5% and phase 
tt/2. 

Data set A contains 100,000 events and spans a period of 3 years, data set B contains 400,000 events and spans 
5 years, finally data set C contains 2,000,000 events and spans 15 years. Events are generated uniformly in 
UTC time and azimuth and with a zenith distribution proportional to sin x cos 8 corresponding to a cosmic 
ray surface array located at longitude 69.25 West and latitude 35.25 South. The events are then accepted using 
an acceptance/rejection method according to the desired combination of modulations. In table 1 a summary of 
the parameters of those various data sets is given. 

On figure [T] the FEMT power as a function of frequency is shown around the diurnal frequency for data set 
A, B and C including the 1% dipolar modulation. For data set A only the diurnal modulation is visible with 
a reduced amplitude of 4% due to the use of the modified time. Neither the sidereal modulation nor the anti- 
sidereal side-band are visible. The resolution, given by the width of the diurnal peak, is however sufficient to 
resolve the diurnal and sidereal frequency. With data set B the situation is clearer, the improved statistics and 
resolution now allow to see the sidereal modulation, the anti-sidereal band is still invisible. Finally with data 
set C one can perfectly resolve the 3 modulations. Note that the anti-sidereal band has an amplitude of about 
0.5% as expected (~ 0.2 x 0.05/2). 

On figure [2] we show the power distribution at the sidereal frequency obtained from 1000 random sets of type 
A, B or C (from left to right). The first row of histograms shows the FEMT with unit weights. Dashed (black) 



3 For a signal p with a uniform random phase and a systematical noise of amplitude b with phase the distribution of the Rice parameter 
PO becomes f(po) <x y/p 2 + b? + 2pb cos 9 / (2pb sin 9) with 9 uniform in [0, 2n]. 
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Figure 1. Data sets A, B and C including the 1% dipolar modulation (see text). Top : Fourier analysis of data set A around 
the diurnal modulation (central violet vertical line) . No signal can be seen at the sidereal (right vertical line) nor at the 
anti-sidereal (left vertical line) frequencies. Middle : Same for data set B, the sidereal modulation is visible. Bottom : Data 
set C, both the sidereal and anti-sidereal modulations are visible. 

curves correspond to the case without RA modulation while solid (red) ones correspond to the case with the 
1.5% RA modulation. One clearly sees on data set C the effect of the 0.5% systematic modulation which 
dominate over the noise and gives a power distribution (solid red) almost uniform between 1 and 2%. 

On the second row of histogram the FEMT is evaluated using event weights. For each event with time U the 
weight Wi is given by : 

W { = [1 +Py COS (OJyU ~ (f)y)]~ l (17) 

where p y and <p y are the power and phase measured at the annual frequency (u) y ) using a unit weight FEMT. 
As expected the improvement is only marginal for data set A but very significant for data sets B and C. 
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Figure 2. Left : Power distribution at the sidereal frequency for data set A with no sidereal modulation (dashed black) and 
with the 1.5% RA modulation (solid red).The curves are gaussian fit to the distributions with the parameters indicated. The 
top histogram corresponds to the Fourier transform on the original data set, the middle histogram includes the correction 
for the annual modulation (see text) and the bottom one for both annual an diurnal. Since the statistical noise is dominant in 
data set A the improvement is marginal after the correction. Middle : Same for data set B. Now the improvement is clearly 
visible between the top histogram and the two lower ones. Right : Same for data set C. Again the improvement is clearly 
visible, the separation between the average noise and the 1.5% modulation increase from top to bottom from 3.7<r to 14<r. 
Note that in this case the red curve of the top histogram which corresponds to the power distribution of a 1.5% modulation 
is no longer a Rice distribution because we are in this case dominated by the side-band noise which has a fixed phase. 

Finally on the last row of histograms of figure [2] we show the power distribution using a weighted FEMT 
correcting for both the annual and diurnal modulations : 

Wi = [(1 + p y cos (ujyU - <j)y) x (l + p d cos (u) dU - <t> <*)]" 1 (18) 

where p y and <f> y are as above and pd and cf>d are the power and phase measured at the diurnal frequency (uj d ). 
Again the improvement is significant for data set C. It is remarkable to note that after this correction all data sets 
show a background power distributions corresponding exactly to the expectation from their limited statistics. 



5. Bi-dimensional application 

The analysis depicted in the previous section dealing with a single declination band is quite appropriate to 
analyze dipole moments (I — 1). However, the power given by equation 11 can be derived, if the statistic is 
sufficient, in several distinct declination bands allowing to reconstruct higher order moments. We just give a 
taste of the method using a sky modulated by a 5% Z = 3, m — 1 YJ m . 

On figure [3] left, we show the reconstructed power for 1 1 bands in declination (from -90 to +20). A total 
of 4 x lO^vents were generated using two different zenith coverage : uniform coverage in zenith (6 in the 
range [0, 90] degrees) and a coverage proportional to cos 9 sin 9 in the range [0, 60] degrees). In both case the 
centre of the detector is at location (69.25° W, 35.25° S). This shows that our analysis results are, as predicted, 
independent of the experimental coverage. On the right figure we show the expected modulation that would be 
seen if a 5% was present on the sky for a statistics of 400,000 events over 5 years. 

The power distribution as a function of declination can be reconstructed independently of the coverage because 
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Figure 3. Left : Reconstructed modulation at the sidereal frequency as a function of declination, for a sky with a 5% Y\ 
modulation and N = 4 x 10 6 . In green uniform coverage in zenith, in red coverage (proportional to cos 9 sin 9). The (blue) 
curve represents the genuine input modulation. As can be seen, although not corrected for, the coverage does not influence 
the reconstructed amplitude evolution with declination. Right a typical plot for 400, 000 events, a 5% Yj 3 modulation and 
a coverage proportional to cos 9 sin 9. 

Eqn 1 1 gives the fraction of events that are modulated within that declination band, independently of the relative 
normalization of one band to another. Hence the result is independent of the eventual variation of the coverage 
with the declination. 



6. Conclusion 



We have constructed a modified Fourier transformation which allows to reconstruct the original right ascension 
modulation of a data set without the need to construct the appropriate coverage map. We also showed that this 
technique allows to obtain, for each order m ^ in right ascension, the dependence in declination over the 
range covered by the observations, independently of the coverage function. 



